%% This is essentially fungi_fluocell_detect, without calling any figures
function [stats, colorvec, fg_bw, fg_bd] = image_preprocess(im, path, name, data, varargin)
%% reverse the image before background subtraction
%     par_name = {'fungi_intensity_th','min_area'};
%     [fungi_intensity_th, min_area] = parse_parameter(par_name, default_value, varargin);
    fungi_intensity_th = data.intensity_th;
    min_area = data.min_area; 
    max_area = 20 * min_area;
    
    %% subtract background on counter staining
    rev_im = 255 * uint8(ones(size(im))) - im;
    file_bg = strcat(path,'output/background_', name,'.mat');
    file_bg_value = strcat(path, 'output/background_value.mat');
    if ~strcmp(data.type, 'Scanscope') || ~exist(file_bg_value, 'file'),
        [bw, ~] = get_background(im, file_bg);
        bg_mean = zeros(3, 1);
        for i = 1:3,
            bg_mean(i) = compute_average_value(rev_im(:,:,i) ,bw);
        end;
        if strcmp(data.type, 'Scanscope'),
            save(file_bg_value, 'bg_mean');
        end;
    else % type is Scancope, and exist file_bg_value
        temp = load(file_bg_value); 
        bg_mean = temp.bg_mean;
    end;
    % Background subtraction followed by converting RGB to Index
    for i = 1:3,
        rev_im(:, :, i) = uint8(rev_im(:, :, i) - bg_mean(i));
    end;
    im_prep = sum(rev_im, 3);

    %% Detect the fungi/objects
    % find the global threshold of the image 
    % The maximum of an automatic threshold and an intensity threshold is used. 
    % Objects less than min_area and larger than max_area are both
    % removed.    
    level = my_graythresh(im_prep);
    fg_bw = im2bw(im_prep, level) & (im_prep > fungi_intensity_th);
    % define a minimal area for fungi and get rid of staining dots
    temp = bwareaopen(fg_bw, min_area); clear fg_bw;
    fg_bw = temp;
    cc = bwconncomp(fg_bw, 8);
    fg_prop = regionprops(cc,'Area', 'PixelIdxList');
    fg_keep = (cat(1, fg_prop.Area) < max_area);
    pixel_index_keep = cat(1, fg_prop(fg_keep).PixelIdxList);
    temp = false(size(fg_bw));
    temp(pixel_index_keep) = true; clear fg_bw;
    fg_bw = temp; clear temp;
    clear temp cc fg_prop fg_keep clear pixel_index_keep;
  
    %% Calculate the characteristics of the detected fungi,
    % such as area, color, aspect ratio, etc.
    % The color is the color of the original image.
    
    % Use the watershed method to separate the bulk detection
    % The average color is the original color. 
    temp = my_WatershedSeg(fg_bw, data.type); clear fg_bw;
    fg_bw = temp; clear temp;
    % get the boundary of the detections
    % Find the accurate detection and get the mean value of each channel
    cc = bwconncomp(fg_bw, 8);
    num_objects = cc.NumObjects;
    avg_color = zeros(num_objects, 3); % for original 3-d image
    %avg_color = zeros(num_objects, 1); % for reversed background subtaction image
%    fg_mask = cell(num_objects, 1);
    for i = 1 : num_objects,
%         % Kathy : This is not efficient, I made the modification below
%         % --- 03/24/2015
%         % It is more efficient to use PixelIdexList to index into the image
%         % matrix
%         im_plain = false(size(im_prep));
%         im_plain(cc.PixelIdxList{i}) = true;
%         fg_mask{i} = uint8(im_plain); % VERY INEFFICIENT
%         for j = 1:size(avg_color,2),
%             %get the mean value on each channel
%             avg_color(i, j) = sum(sum(fg_mask{i} .* im(:,:,j))) / length(cc.PixelIdxList{i}); %based on original image           
%         end
%
        pixel_index = cc.PixelIdxList{i};
        %get the mean value on each channel
        for j = 1:size(avg_color, 2),
            temp = im(:, : , j);
            avg_color(i, j) = sum(temp(pixel_index))/length(pixel_index); 
            clear temp;
        end;
        clear pixel_index;
    end; 
    % Get the properties of each detection
    fg_prop = regionprops(cc, im_prep, 'PixelIdxList','Eccentricity','Solidity','Area', 'Centroid');
    % Boundary roughness
    fg_bd = bwboundaries(fg_bw, 'noholes');
    fg_bd_roughness = zeros(size(fg_bw));
    boundary_roughness_vector = zeros(num_objects, 1);
    for i = 1:num_objects,
        % first direvative along the x and y, since fb_bd{i} is nx2
        boundary_roughness_vector(i) = sum(sum(abs(diff(fg_bd{i}, 1, 1))))/...
            sqrt(fg_prop(i).Area);
        fg_bd_roughness(cc.PixelIdxList{i}) = boundary_roughness_vector(i);
    end;

    % the following assumes that there are actually detections
    if(not(isempty(fg_prop)))
        fg_keep = (cat(1, fg_prop.Eccentricity) <= data.eccentricity_th) & ...
            (cat(1, fg_prop.Solidity) >= data.solidity_th)& ... 
            boundary_roughness_vector <= data.roughness_th;
        pixel_index_keep = cat(1, fg_prop(fg_keep).PixelIdxList);
        temp = false(size(fg_bw));
        temp(pixel_index_keep) = true; clear fg_bw;
        fg_bw = temp; clear temp;
        clear fg_prop fg_bd num_objects;
        stats = regionprops(fg_bw, im_prep, 'Area', 'Eccentricity', 'MeanIntensity',...
            'MajorAxisLength', 'BoundingBox', 'Centroid');
        fg_bd = bwboundaries(fg_bw, 'noholes');
    else % in the case that there are no detections
        clear fg_prop fg_bd num_objects;
        stats = regionprops(0, 'all'); % we still want this to have the regionprops struct
        fg_bd = bwboundaries(fg_bw, 'noholes');
    end
    colorvec = avg_color;
return;

